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Abstract 

We develop numerical methods to calculate electron dynamics in crystalline solids in real-time 
time-dependent density functional theory employing exchange-correlation potentials which repro¬ 
duce band gap energies of dielectrics; a meta generalized gradient approximation (meta-GGA) 
proposed by Tran and Blaha [Phys. Rev. Lett. 102, 226401 (2009)] (TBm-BJ) and a hybrid 
functional proposed by Heyd, Scuseria, and Ernzerhof [J. Ghem. Phys. 118, 8207 (2003)] (HSE). 
In time evolution calculations employing the TB-mBJ potential, we have found it necessary to 
adopt a predictor-corrector step for stable time-evolution. Since energy functional is not known 
for the TB-mBJ potential, we propose a method to evaluate electronic excitation energy without 
referring to the energy functional. Calculations using the HSE hybrid functional is computation¬ 
ally expensive due to the nonlocal Fock-like term. We develop a computational method for the 
operation of the Fock-like term in Fourier space, for which we employ massively parallel computers 
equipped with graphic processing units. To demonstrate significances of utilizing potentials pro¬ 
viding correct band gap energies, we compare electronic excitations induced by femtosecond laser 
pulses using the TB-mBJ, HSE, and a simple local density approximation (LDA). At low laser 
intensities, electronic excitations are found to be sensitive to the band gap energy: results using 
TB-mBJ and HSE are close to each other, while the excitation of the LDA calculation is more 
intensive than the others. At high laser intensities close to a damage threshold, we have found that 
electronic excitation energies are similar among the three cases. 
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I. INTRODUCTION 


At current frontiers of optical sciences, interactions of intense and ultra-short laser pulses 


with solids are attracting much interests [l-3|. For example, ultra-short laser pulses with 
intensities above a damage threshold are used for non-thermal laser-processing [4,1^. Ultra- 
short laser pulses irradiated on dielectrics close to the damage threshold have been found to 
show iutnguiug phenomena such as producng aa ultrafast charge trarrsfer g and modifying 
the band gap in sub-femtosecond time scale [71] . These phenomena induced by ultra-short, 
intense laser pulses are expected to open up new technological applications. 

For microscopic understanding of these phenomena, it is quite significant to theoretically 
investigate nonlinear electron dynamics in the medium excited by the strong electric field of 
laser pulses. Real-time calculation based on the time-dependent density functional theory 


(TDDFT) |8(| is a useful theory for the investigation of nonlinear electron dynamics. It has 
been applied to nonlinear laser-solid interactions such as optical breakdown 9|, generation 


of coherent phonon 


10| . and calculation of nonlinear optical properties [n|, as well as linear 


responses for a weak field 121. 

In order to apply the TDDFT to practical problems, we need to use approximate 
exchange-correlation potentials. In most applications of the TDDFT, adiabatic approx¬ 
imation, in which the exchange-correlation potential at each time is determined by the 
density (, current, and kinetic density if necessary) at the same time, is assumed. In 
the following developments, we will consider three exchange-correlation potentials: local 
density approximation (LDA), meta generalized gradient approximation (meta-GGA), and 
hybrid functional, and use them in both ground state and time evolution calculations in the 
adiabatic approximation. 

Among the three approximations, the LDA is the simplest one. We will use the energy 
functional of Ref. 13|. Because of the well-known band gap problem in the static LDA 
calculations, optical gaps of semiconductors and insulators are also underestimated in the 
TDDFT using the adiabatic LDA. The underestimation of the optical gap causes serious 
problems when we compare calculated physical quantities with experimental results in non¬ 
linear electronic excitations induced by intense laser pulses as well as in linear responses for 
a weak optical fields. 

Recently, Tran and Blaha proposed a meta-GGA exchange potential which systematically 


2 














improves the band gap energy of various materials 


Q- 


We abbreviate the potential as 


TB-mBJ potential. Calculated band gap energies using the TB-mBJ potential agree with 
experimental gap energies in a quality similar to those by the GW method jl^ . 


Hybrid functional is, at present, one of the most successful approximations in practical 
applications of static and time-dependent density functional theory. In solid state physics, 
a hybrid functional that is proposed by the Heyd, Scuseria, and Ernzerhof (HSE) has been 
widely applied 1^ . The HSE functional signihcantly improves fundamental gap comparing 
to one in the LDA. It has been shown that also accurate optical absorption spectra in several 
dielectrics and partial inclusion of excitonic excitation by J. Paier et al 16|. 


In this work, we present numerical methods to achieve time-evolution calculations in the 
real-time TDDFT employing the TB-mBJ and HSE potentials. We have already reported re¬ 


sults employing the TB-mBJ potential in t 


in both linear and nonlinear regimes 


17 


le real-time TDDFT for light-matter interactions 
19| . To carry out stable time-evolution calculations. 


we hnd that a predictor-corrector step is indispensable in the calculations using the TB-mBJ 
potential, while it is not necessary in the calculations using the LDA and HSE potentials. 
The TB-mBJ potential is given in Ref. without referring to the energy functional, and 
an explicit form of the energy functional providing the TB-mBJ potential is not known. 
We discuss how to evaluate electronic excitation energies during the irradiation of a laser 
pulse without referring to the energy functional. In order to calculate nonlinear electron 
dynamics using the HSE functional, we develop a Fourier-space method for the Fock-like 
operations: the non-local Fock-like term in the HSE functional is calculated in Fourier space, 
making use of an efficient Fast-Fourier-Transformation library on an accelerator type super¬ 
computer. We apply the above methods for crystalline silicon under irradiation of ultrashort 
laser pulses, and compare electronic excitations by the three potentials. 


The paper is organized as follows. In Sec. H, we describe methods of electron dynamics 
calculations in crystalline solids utilizing the TB-mBJ and HSE potentials. In Sec. HI and 
IV, we present calculations for electron dynamics in crystalline silicon, using LDA, TB-mBJ, 
and HSE potentials in linear and nonlinear regimes, respectively. In Sec. V, a summary is 
presented. 
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II. FORMALISM 


A. Electron dynamics calcnlation in crystalline solids 


First we briefly explain a general feature of electron dynamics calculations in crystalline 
solids based on the TDDFT. Details of the method are explained in Ref. 12|. We consider 


electron dynamics in crystalline solids under visible or infrared laser pulses. A typical 
wavelength of the electric held of such laser pulses is a micrometer, while the length scale of 
the electron dynamics induced by the laser pulse is less than a nanometer. Therefore, in a 
unit cell of crystalline solids, we can treat the laser electric held as a spatially-uniform and 
time-dependent electric held E{t). This is a dipole approximation for interactions between 
the laser pulse and electrons in a medium. The electron dynamics is described by the 
following time-dependent Kohn-Sham (TDKS) equation, 

( 1 ) 

where hxs{t) is the time-dependent Kohn-Sham Hamiltonian. It is given by 


d 


hRsit) — 


2mp 


^ion 


( 2 ) 


where A{t) is the vector potential which is related to the electric held E{t) by A{t) = 
—cf*dt'E{t'). The Hamiltonian contains the Hartree potential f), an exchange- 

correlation potential nxc(Df), and an ionic potential Vion[A{t)]. The Hamiltonian of Eq. ([2]) 
is periodic in space so that we may apply the Bloch theorem to the orbitals 'ipiifif) each 


time. For the ionic potential, we employ a norm-conserving pseudopotential 


20 1 with the 


separable approximation [2l|. The ionic potential depends on the vector potential A{t) due 


to the nonlocality of the potential 


m 


The electric held E{t) may include an induced polarization held which depends on macro¬ 
scopic geometry of the material which we treat. In the present calculation, we assume the 


transverse geometry explained in Ref. 


22l | in which the induced polarization helds do not 


appear. Moreover, the vector potential A{t) may include exchange-correlation terms in the 


time-dependent current density functional theory 23 


24| |. We simply ignore them, since 


reliable potentials have not yet been available at present. 
We dehne a matter current density of electrons j(F, t). 


J(f, t) = — ^ 3f? 


nip 


iJ*i{ryt){-inV + -A{t)}il)i{fyt) 
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( 3 ) 
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and spatially averaged electric current density J(t)., 

At) = -e^ / j{f,t)df+ JNL{t), 
iZ Jn 


( 4 ) 


where hZ is a volume of the unit cell and JNbit) is a contribution from the nonlocal part of 
the ionic pseudopotential |^. 


B. Exchange-correlation potentials 

We will use the LDA, TB-mBJ, and HSE potentials in the real-time TDDFT calculations 
in the adiabatic approximation. We present below brief explanations for the TB-mBJ and 
HSE potentials. 

The TB-mBJ exchange potential in the adiabatic approximation has the following form, 


VTB-mBj{r,t) = c^VBR{r,t) + (3c^ - 2)- 


_ rjf, t) 

TT V 12^ 


( 5 ) 


where VBRif^t) is the Becke-Roussel potential |25j, p{r,t) = is the electron 

density, and T{r,t) is the kinetic energy density given by 


r r 






{-iW -A{t)}ipi{f,t) 
c 


( 6 ) 


It has been known that the band gap depends crucially on the parameter Cm of Eq. 


26l |. In the original TB-mBJ paper [1^, a formula for the mixing parameter is proposed 


using the electron density and its gradient. In this work, we treat the parameter Cm as 
adjustable and determine the value in such a way that the band gap coincides with the 
experimental gap of silicon. As will be shown, we may hnd the value of Cm that reproduces 
the experimental optical gap as well as the band gap simultaneously. The potential of Eq. 
(jS]) is not gauge invariant, since the kinetic energy density of Eq. (|6]) is not 


27 


30|. To 


recover the gauge invariance, we replace all kinetic energy density terms in the potential as 

r(r, t) r(r, t) - j^{r, t)/p{f, t), (7) 


where j{r,t) is the current density of Eq. 

A91 correlation potential M|| following t 
The HSE exchange-correlation potential is given by 


]). For the correlation potential, we use the 

n 

PW91 correlation potential [M[ following the original TB-mBJ paper |14| . 


_ -- _ }i.,SR,PBE , }^ ^SR,Fock 

^xc — ^xc ^ 4 ^ ’ 


PBE 


( 8 ) 
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where is the PBE exchange-correlation potential 32|, and ^.pg 

PBE exchange and the non-local Fock-like potentials corresponding to the short-range part 
of the Coulomb potential. The HSE exchange-correlation potential includes a parameter /i, 
the range separation parameter for the Coulomb interaction. In this work, we use the value 
fi of 0.3 A“h 


C. Numerical method 


We solve the TDKS equation, Eq. ([T]), in real time and real space |33|, [3^. We choose a 


cubic unit cell containing eight silicon atoms and 32 valence electrons in total. We consider 
electron dynamics only, freezing atomic positions at their equilibrium positions in the ground 
state. The unit cell is discretized by 16^ uniform grid points. The first Brillouin zone is 
discretized by 24^ fc-points for calculations using the LDA and TB-mBJ potentials, and 12^ 
fc-points when using the HSE potential. We use the preudopotential which was constructed 
using the LDA. 

1. Predictor-corrector procedure 


For time-evolution calculations, we repeat a time evolution of a small time step At: 

At\ 


ipiir, t At) ^ exp 


hKs{t + 2 
ih 


-At 




(9) 


where the Hamiltonian is approximately evaluated at a time t At/2. We expand the time 


33|,l35|: 


evolution operator, exp ^)At/ifi , in Taylor series and take up to fourth order 


exp 


At' 


hKsjt + 2 
ih 


-At 


^ 1 

F- 

^ rt' 

n=0 


hKs{t -\- 


At' 


ih 


-At 


( 10 ) 


In the above procedure, we need to guess the Hamiltonian at a time t At/2. In calcu¬ 
lations using the LDA and HSE potentials, we find that use of the Hamiltonian at time t, 
hxsii + At/2) ~ hxsit), provides results with sufficient accuracy. However, we have found 
that the predictor-corrector procedure to approximately evaluate hxsit + At/2) is indis¬ 
pensable for stable time-evolution when we employ the TB-mBJ potential. The predictor- 
corrector procedure is usually carried out in the following steps 


We first calculate wave 
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( 11 ) 


functions t + At) using the Hamiltonian at time t, 


i + At) = exp 


hKs{t)^t 

ih 


'll^i{f,t). 


Using the Kohn-Sham orbitals we next construct the Hartree potential and the 

exchange correlation potential at time t + At. We then approximate the Hamiltonian 
at time t + At/2 by averaging the potentials at time t and t + At as follows: 


^ ’ 2 2 


. I t + At) + vUa t) 

^XC\’ 5 t, -f* 2 / ^ 2 


( 12 ) 

(13) 


Finally, we calculate the Kohn-Sham orbitals at t-|-At by Eq. (5) using the Hamiltonian at 
time t -|- At/2. 

In the following, we numerically demonstrate the necessity of the predictor-corrector 
procedure when the TB-mBJ potential is used. Figure [U shows the electric current in the 
crystalline silicon as a function of time after an instantaneous weak distortion is applied 
at t = 0. A time step of At = 0.04 a.u. is used. Figure [U (a) shows the electric current 
using the LDA potential, and (b) shows the electric current using the TB-mBJ potential. In 
both panels, red-solid lines show results using the predictor-corrector procedure, while green- 
dashed lines show results without using the predictor-corrector procedure. One sees that, 
while the calculations using the LDA potential give the same result regardless of whether the 
predictor-corrector procedure is used, the calculations using the TB-mBJ potential proceed 
stably only if the predictor-corrector step is used. Until t = 0.8 fs (about 800 time steps), the 
calculated electric currents with and without the predictor-corrector procedure coincide with 
each other. However, at t = 0.8 fs, the calculation without the predictor-corrector procedure 
starts to show unphysical oscillations. We carried out calculations with and without the 
predictor-corrector procedure for several time steps. At = 0.02, 0.04, and 0.08 a.u., and have 
found that results without the predictor-corrector procedure always fail showing unphysical 
oscillations after a certain number of iterations. Calculations using the predictor-corrector 
procedure give the same result for the three time steps. We thus conclude that the predictor- 
corrector procedure is indispensable for stable time-evolution of orbitals when using the TB- 
mBJ potential. In the calculations presented below, we use the predictor-corrector procedure 
when we employ the TB-mBJ potential. 
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FIG. 1. (color online) Electric currents as functions of time after an instantaneous weak distortion 
is applied at t = 0. Panel (a) shows the results using the LDA potential, and panel (b) shows 
the results using the TB-mBJ potential. In both panels, calculated results with (red-solid) and 
without (green-dashed) the predictor-corrector procedure are compared. 

2. Nonlocal Fock-like term 

In calculations using the HSE potential, we need to calculate operations of nonlocal 
Fock-like term on orbitals. The nonlocal Fock-like term y^NFcok following form, 

j=occ 

where only occupied orbitals are summed. U 5 _R(r) is the short-ranged Coulomb-like potential, 
= erf(yU,r)/r. In practical calculations, we treat the nonlocal Fock-like term in Fourier 
space as: 

E I dGe-“^-<’v‘”‘(G)p,i(G,t), 

j=OCC ^ ^ ' 

where v^^{G) and pji{G,t) are Fourier transformation of and orbital products 

Pji{f,t) = 'ijj*{r,t)'ijji{r,t), respectively. 

The Fourier space operation has been carried out efficiently employing the Fast-Fourier- 
Transformation library on large-scale accelerator type supercomputers. In the HSF calcula¬ 
tion, we further employ so-called down-sampling method 0 to reduce the computational 
cost, which we briefly describe below. First, we decompose the orbital index i into the 
reciprocal lattice vectors of the unit cell k and the band index b: Using 

this notation, we rewrite Fq. flTbll as follows: 

= Y. j dGe-''^^^v^^{G)p^,^^^^{G,t), (16) 

q^b'=occ ^ ' 

























where the orbital product is described as P^i using the new in¬ 

dices. In the down-sampling method, the summation of Eq. ffT6|) is carried out reducing the 
number of (f-points. In this work, we sample 12^ fc-points to describe Kohn-Sham orbitals, 
while we sample 4^ g-points to evaluate the operation of the non-local Fock-like term of Eq. 

m- 


III. LINEAR RESPONSE CALCULATION IN REAL TIME 


In this and next sections, we apply the real-time TDDFT calculations for electron dy¬ 
namics in crystalline silicon using the LDA, TB-mBJ, and HSE potentials. Although linear 
response properties of solids may be investigated either in frequency or in time domains, 
time-domain calculation is a unique option for strongly nonlinear dynamics induced by ultra- 
short, highly-intense laser pulses. In this section, we discuss a real-time calculation for a 
dielectric function of silicon. The purpose of this section is to demonstrate that numeri¬ 
cal methods of time evolution calculations presented in the previous section work stably in 
real-time linear response calculations. 

For dielectric functions, several ab-initio approaches have been developed. Among them, 
GW plus Bethe-Salpeter equation approach |38|, l39|, which is based on many-body per¬ 
turbation theory, has been successful to describe dielectric functions of various materials 
quantitatively. Ah-initio calculations based on TDDFT have also been presented [^. Re¬ 
cently, dielectric functoins in the TDDFT have been reported using meta-GGA H and HSE 
functionals, which are successful in describing band gap energies of various materials 
correctly. 

We briefly explain how we calculate dielectric function as a function of frequency in real¬ 
time TDDFT calculations. We calculate time evolution of electron orbitals under a weak 
electric field, E{t) = —solving Eq. ([1]). Using electric conductivity crij{t), the 
spatially-averaged electric current J(t), which is given by Eqs. ([3]) and (j4]), is related to the 
electric held E(t) by 


Mt) J - t')Ej{t'), (17) 

j 

where i and j indicate Gartesian component, i,j = x,y,z. Taking Fourier transformations 
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of both sides, we obtain a relation in frequency representation, 


Ji{^) = Y. ( 18 ) 

j 

where Ei{u)), Ji{oj), and (Tij^oj) are Fourier transformations of Ei(t), Ji(t), and Cijit), re¬ 
spectively. Frequency-dependent dielectric function eij{u) is related to as usual, 

47r2 

UJ 

In principle, we may use any electronic field Ei{t) to calculate For numerical 

convenience, we employ a 5-type distortion in 2 ;-direction for the electric field, E{t) = sCzSit) 
where s is a small parameter. Then the conductivity is proportional to the electric current, 
o'zz{t) = Jz{t) / s- This relation tells us that the Fourier transformation of the electric current 
provides the frequency-dependent conductivity. To reduce numerical noises caused by the 
finite evolution time T, we use a mask function M{x) in the Fourier transformation as 
follows. 


(TzzAOJ] = 


s Jo 


dtJz{t)E^^M{t/T). 


( 20 ) 


In practice, we employ the mask function M{x) = 1 — 3x^ + 2x^ 421. 

Figure [2] shows electric currents Jz(t) in crystalline silicon as functions of time after the 
delta-function distortion is applied at t = 0. Red-solid line shows the result using the HSE 
potential, green-dashed line shows the result using the TB-mBJ potential, and blue-dotted 
line shows the result using the LDA potential. All currents accurately coincide with each 
other at t = 0, since the initial current is determined by the sum rule |^. Shortly after 
the distortion (0 fs to 1 fs), all currents look similar to each other. They start to depart 
after a few femtoseconds. As discussed above, the electric current is proportional to the 
conductivity as a function of time when we employ the delta-type distortion. Therefore, the 
difference of the currents directly reflects the difference of the conductivities. 

From the currents shown in Fig. [2], dielectric functions can be calculated using Eqs. (IT^ 
and fl20l) . Figure [3] shows the results: the real part in (a) and the imaginary part in (b). 
Red-solid line shows the result using the HSE potential, green-dashed line shows the result 
using the TB-mBJ potential, and blue-dotted line shows the result using the LDA potential. 
An experimental result 43(] is also shown as black dash-dotted line. Our results using the 
HSE and TB-mBJ potentials shown in Fig. [3] almost coincide with previous calculations 


using frequency-domain methods 


16 


4l|. 
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FIG. 2. (color online) Electric currents in crystalline silicon as functions of time after an impulsive 
distortion is applied at time t = 0. Results using the three exchange-correlation potentials are 
compared: HSE (red-solid), TB-mBJ (green-dashed), and LDA (blue-dotted). 

As seen from the panel (b), the optical gap is underestimated when the LDA potential is 
used. This is closely related to the well-known band gap problem of the LDA. In contrast, 
the optical gap is well reproduced when we use the HSE and TB-mBJ potentials. The 
real parts of the dielectric functions calculated using the HSE and TB-mBJ potentials also 
reproduce accurately the experimental values in low frequency region (below 3 eV). 


IV. NONLINEAR ELECTRON DYNAMICS UNDER INTENSE LASER PULSES 

In this section, we investigate nonlinear electron dynamics in crystalline silicon under 
intense, ultra-short laser pulses using the three exchange-correlation potentials. We employ 
electric helds derived from the following vector potential, 

[ — ^ cos(a;t) sin^(^t) (0 < t < T^) 

A{t) = J ^ ^ ^ V dj ^21) 

I 0 (otherwise), 

where ca is a mean frequency of the held, is the pulse duration, and Eq is the maximum 
electric held of the pulse. We set u and to 1.35 eV/h and 16 fs, respectively. We dehne 
the maximum intensity of the laser pulse Jq as Jq = cEq/Stt. We calculate electron dynamics 
using three exchange-correlation potentials changing the intensity Jq. We will examine two 
physical quantities: the excitation energy during and after the laser irradiation, and the 
number density of excited electrons after the laser irradiation. 
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FIG. 3. Dielectric functions of silicon; real part in (a) and imaginary part in (b). Results using 


the three exchange-correlation potentials are 
and LDA (blue-dotted). An experimental result 


compared: HSE (red-solid), TB-mBJ (green-dashed), 
43l ] is also shown as black-dash-dotted lines. 


A. Excitation energy 


We first consider electronic excitation energy induced by laser irradiation. The excitation 
energy is one of the most important quantities to investigate laser-matter interactions. For 
example, optical damage thresholds can be estimated by comparing the electronic excitation 


energy with cohesive energy of solids 


M, il- 


Although the LDA and HSE potentials are 


explicitly derived from energy functionals, the TB-mBJ potential is not derived from any 
functionals but given directly. It is even unclear whether there exists an energy functional 
which provides the TB-mBJ potential. We first discuss how to evaluate electronic excitation 
energy when we use the TB-mBJ potential. We then investigate differences in electronic 
excitation energies using the three exchange-correlation potentials. 

We first note that there are two procedures to calculate electronic excitation energy which 
should give the same results if the exchange-correlation potential is derived from an energy 
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functional. One is to use the explicit expression of the energy functional given by 


Eex{t) =J 2 


(-ifiV + lAit)) 

2m 


+ Vion[A{t)] 


Mr,t) 


+ ^I dfvnir, t)p{r, t) + t)}] - Eqs, 


(22) 


where E^c is the exchange-correlation energy functional and Eqs is the energy of the gronnd 
state. The other is to calculate the work done on electrons by the laser electric held, 

W{t)= dt'J{t') ■ E{t'), (23) 

J —oo 

where J{t) is the electric current dehned by Eq. (|1]). When the exchange-correlation po¬ 
tential is derived from the energy functional Eex[{'^i{T,t)}]j we may easily prove that there 
holds 

y^it) = f dt'^E,^{t') = E,^{t). (24) 

J —OO CLt 

In calculations using the LDA and HSE potentials, we may employ either Eq. fl2^ or Eq. 
(|23|) . In calculations using the TB-mBJ potential, Eq. (|23|) is the unique option. 

Althongh Eq. (12^ and Eq. (|23ll are analytically equivalent to each other if the potential 
is derived from an energy functional, we have found that Eq. fl22ll is numerically more 
favorable when the applied laser pulse is weak. Figure |4] shows E^xif) of Eq. fl22D using the 
LDA potential as a function of time for two intensities, I = 1.0 x 10^*^ W/cm^ in (a) and 
J = 1.0 X 10^^ W/cm^ in (b). In the panel (a), we hnd an appreciable excitation energy 
during irradiation of the laser pulse. However, at the end of the pulse, the system almost 
returns to the ground state since the laser frequency is lower than the direct band gap. 
We shall call it virtual excitation. If we calculate the excitation energy using Eq. (12^ . 
the integrand J{t) ■ E{t) behaves oscillatory during the irradiation of the laser pulse. The 
integration over time mostly cancels and gives the total work done by the laser pulse which 
is very small as seen from the panel (a). Therefore, the calculation using Eq. fl23|) may suffer 
from an accnmulation of a numerical error during the time integral. Figure 0] (c) compares 
the electronic excitation energies after the weak laser pulse ends (Jq = 1.0 x 10^°W/cm^) 
using several different time steps At. Red up-pointing and green down-pointing triangles 
show the excitation energies calculated by Eq. fl22D and by Eq. ([23]), respectively. As seen 
from the hgure, results coincide numerically if we employ sufficiently small time step. The 
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Time (fs) At (a.u.) 

FIG. 4. Laser-induced electronic excitation energy calculated using the LDA potential. Panels 
(a) and (b) show the electronic excitation energy (red-solid) and the applied squared electric field 
(green-dashed) as functions of time. Panels (c) and (d) show the electronic excitation energy after 
the laser irradiation ends using Eqs. (I22|) and (1231) . Results in the panels (a) and (c) are for the 
irradiation of the weak laser pulse (Jq = 1.0 x 10^'^W/cm^), while results in the panels (b) and (d) 
are for the irradiation of the strong laser pulse (Iq = 1-0 x lO^^W/cm^). 

error in the excitation energy increases linearly with the time step At when we employ Eq. 

When the laser pulse is so strong that electrons are substantially excited by the laser 
pulse, we hnd a steady increase of the excitation energy in time, as seen in Fig. H] (b). We 
shall call it real excitation. In this case, the difference in the numerical calculation is very 
small between the calculations using Eqs. (12^ and (|23ll . even if we employ the largest time 
step. At = 0.08 a.u. 

When we use the TB-mBJ potential, we can only evaluate excitation energy using Eq. 
([23]). To obtain excitation energy, we carry out two calculations using different time steps. 
At = 0.02 a.u. and At = 0.01 a.u. We then estimate the converged excitation energy by a 
linear extrapolation using the two calculations. 
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FIG. 5. (color online) Electronic excitation energies as functions of time induced by laser irradia¬ 
tion. Panel (a) shows the excitation energy under irradiation of a weak laser pulse (Jq = 1.0 x 10^^ 
W/cm^), while panel (b) shows that of a strong laser pulse (Iq = 5.0 x 10^^ W/cm^). Results using 
the three exchange-correlation potentials are compared: HSE (red-solid), TB-mBJ (green-dashed), 
and LDA (blue-dotted). Squared electric fields are also shown as gray-dash-dotted lines. 

Figure E] shows electronic excitation energies as functions of time using the three exchange- 
correlation potentials; LDA, TB-mBJ, and HSE. Figure |5] (a) shows results for a weak laser 
pulse (Jo = 1.0 X lO^^W/cm^), and Figure [5] (b) shows results for a strong laser pulse 
(Jo = 5.0 X lO^^W/cm^). Squared electric fields E‘^{t) are also shown as gray dash-dot 
lines in each panel. As seen from Fig. [5] (a), the electric field induces virtual excitation 
during the laser irradiation. After the irradiation ends, the excitation energies using the 
HSE and TB-mBJ potentials return to almost zero. In contrast, the excitation energy using 
the LDA potential remains finite after the irradiation, showing a substantial real excitation. 
These results indicate that wider optical gap in the calculations using the TB-mBJ and HSE 
potentials suppresses the real excitation (see Fig. ED. 

Figure [5] (b) shows that electronic excitations induced by the strong electric field persist 
even after the pulse ends irrespective of the exchange-correlation potentials. One sees that 
the excitation energies using the HSE and TB-mBJ potentials are similar to each other, while 
the excitation energy using the LDA potential is larger than the others at each time. This 
is consistent with the above finding that the wider band gap suppresses the real excitation. 

Next we investigate the electronic excitation energy after laser pulses end, changing the 
laser intensity Jq. Figure [6] shows excitation energies after the pulses end {t = 19.3 fs). Red 
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FIG. 6. (color online) Excitation energies as functions of laser intensity calculated using the three 
exchange-correlation potentials: HSE (red down-pointing triangle), TB-mBJ (green up-pointing 
triangle), and LDA (blue square). Cubed and squared intensities normalized to the excitation 
energies at low intensity are also shown as black-solid and gray-dashed lines, respectively. 


down-pointing triangle shows the HSE resnlt, green np-pointing triangle shows the TB-mBJ 
resnlt, and bine sqnare shows the LDA result. Lines of cubed and squared intensity of the 
laser pulses are also shown as black-solid and gray-dashed lines, respectively. These lines 
are normalized at low intensity to coincide with the excitation energies. Black horizontal 
line shows the cohesive energy of silicon, 4.62eV/atom which may be regarded as a 
reference of the damage threshold. 


Multi-photon excitations are expected for weaker intensities, while tunneling excitations 


are expected for stronger intensities ^ . At low intensity region, the excitation energy using 
the LDA potential can be £t with the squared intensity line Jq, while the excitation energies 
calculated using the TB-mBJ and HSE potentials can be £t with the cubed intensity line 
Jq. These behaviors indicate that two and three-photon absorption processes take place 
at the low intensity region for the LDA and the others cases, respectively. The numbers 
of absorbed photons are consistent with the ratio of the photon energy of the laser pulse 
{Tioj =1.35 eV) to the optical gap: the optical gap of silicon is 2.5 eV in the LDA case, and 
it is 3.1 eV in the HSE and the TB-mBJ cases (see Fig. [3]). 


One sees that the excitation energy using the LDA potential is much higher than the 
excitation energies using the HSE and TB-mBJ potentials at low intensity region. In con¬ 
trast, the difference among the excitation energies becomes relatively small at high intensity 
region close to the cohesive energy of the medium. 
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B. Number density of excited electrons 


The number density of excited electrons due to appliec 


observables to characterize laser excitation processes 


48 


ectric fields is one of significant 


49| . To define the number density 


of excited electrons after laser irradiation, we first define eigenstates of the Kohn-Sham 
Hamiltonian after the laser pulse ends (t/ = 19.3fs), 




(25) 


Using these eigenstates, we define the number density of excited electrons riex by, 

nex = nelec-^ I (^/)) T) ( 26 ) 

^ i,j=occupied 

where Ueiec is the number density of valence electrons in the ground state, hi is the volume of 
the unit cell, and 'ipi{f,t) are the time-dependent Kohn-Sham orbitals, which are solutions 
of Eq. ([1]). Only occupied states are summed in Eq. fl26|) . 

Figure [7] shows the number density of excited electrons as functions of laser intensity Jq. 
We find similar features to those seen in excitation energies shown in Fig. |6l The number 
density using the HSE potential is very close to the number density using the TB-mBJ 
potential. The number density using the LDA potential is larger than the others. While 
the number density using the LDA potential is proportional to the squared intensity of laser 
pulse in the low intensity region, the number densities using the HSE and TB-mBJ potentials 
are proportional to the cubed intensity. At high intensity region, the number densities of 
excited electrons using three exchange-correlation potential are similar to each other. 

Figure [H] shows excitation energy per excited electron, namely, the excitation energy 
shown in Fig. |6] divided by the number of excited electrons shown in Fig. [7] as a function 
of laser intensity Jq. Twice and three times of the photon energy, 1.35 eV, are also shown 
by horizontal lines. At low intensity region, the excitation energy per excited electron using 
the TB-mBJ and HSE potentials approach to three times of the photon energy, 3 x 1.35 eV, 
while that using the LDA potential approaches twice of the photon energy, 2 x 1.35 eV. This 
fact obviously indicates that two and three-photon absorption processes dominate in the 
LDA case and the other cases, respectively, and is consistent with the intensity dependences 
seen in Figs. |6]and[71 
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FIG. 7. (color online) The number densities of excited electrons as functions of laser intensity calcu¬ 
lated using the three exchange-correlation potentials: HSE (red down-pointing triangle), TB-mBJ 
(green up-pointing triangle), and LDA (blue square). Cubed and squared intensities normalized 
to the number densities of excited electrons at low intensity are also shown as black-solid and 
gray-dashed lines, respectively. 
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FIG. 8. (color online) Excitation energy per excited electron is plotted against the laser intensity 
calculated using the three exchange-correlation potentials: HSE (red down-pointing triangle), TB- 
mBJ (green up-pointing triangle), and LDA (blue square). Two horizontal lines are shown at 
2 X 1.35 (black-dashed) and 3 x 1.35 eV (gray-dotted). 

V. SUMMARY 

In this paper, we developed methods to carry out real-time TDDFT calculations using the 
TB-mBJ meta-GGA potential and the HSE hybrid functional, which are known to reason¬ 
ably reproduce band gaps of insulators. To carry out stable time evolution of orbitals using 
the TB-mBJ potential, we found that the predictor-corrector procedure is indispensable. 
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Without the predictor-corrector procedure, the electric current induced by an impulsive ex¬ 
citation shows unphysical oscillations after the time evolution of a few femtoseconds. This 
failure cannot be avoided even if we employ very small time step. In the calculations using 
hybrid functional, we developed a Fourier-space method to operate Fock-like term of hy¬ 
brid functional on Kohn-Sham orbitals, making use of recently developed accelerator-type 
supercomputers. 

We examined electron dynamics using the TB-mBJ, HSE, and LDA potentials. In linear 
response calculations, we investigated the dielectric response of crystalline silicon in time- 
domain. Calculated dielectric functions using the HSE and TB-mBJ potentials are conhrmed 
to be similar to the previous results of frequency-domain TDDFT calculations 
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4ll. In 


electron dynamics calculations under strong, untra-short laser pulses, we investigated the 
excitation energy and the number density of excited electrons during and after the laser 
irradiations. Since an explicit form of the energy functional corresponding to the TB-mBJ 
potential is not known, we developed a method to evaluate electronic excitation energy, 
not referring to the energy functional but calculating a work done by the external held to 
electrons. For irradiation of weak laser pulses, results using the HSE and TB-mBJ potentials 
are almost the same to each other, while the excitation energy and the number density of 
excited electrons of the LDA calculation are much larger than the others. For irradiation of 
strong laser pulses close to the damage threshold, we found that results are similar among 
the three potentials. 

We consider that real-time TDDFT calculations using sophisticated exchange-correlation 
potentials such as the TB-mBJ and HSE potentials, which describe the band gap of various 
insulators reasonably, will enables us to carry out quantitative descriptions for electron 
dynamics in crystalline solids even under extremely nonlinear conditions. They are expected 
to provide signihcant insights and knowledge in wide helds of optical sciences including 
nonlinear optical responses, optical-control of electrons, and laser processing in solids. 
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